%% Michael D. Konig, Andrei Levchenko, Tim Rogers and Fabrizio Zilibotti
%% Aggregate fluctuations in adaptive production networks
%% Forthcoming in Proceedings of the National Academy of Sciences, 2022

%% Replication file (Matlab) for Figure 5, main text.

%% Requirements: 

% "firm_factset_ids_sic_nation.csv" -> CSV file with the following colums: id; firm_factset_id; firm_sic; firm_nation
% "2015_adjacency_matrix.mat" -> Adjacency matrix (array) stored in Matlab data file (for the year 2015).
% "mle_eta.mat" -> Matlab file with estimated firm fixed effects. Outpt generated from "mle_attachment.m".
% "in_deg_years_us_manuf.mat" -> Matlab data file with firms' in-degrees/number of suppliers (US manufacturing firms only).
% "in_deg_years_us.mat" -> Matlab data file with firms' in-degrees/number of suppliers (US firms only).
% "in_deg_years_manuf.mat" -> Matlab data file with firms' in-degrees/number of suppliers (manufacturing firms only).
% "in_deg_years.mat" -> Matlab data file with firms' in-degrees/number of suppliers (full sample).
% "European_countries.csv" -> CSV file with a list of European countriesand the following colums: Name,2 Letter Country Code,3 Letter Country Code,EU

clear all;
close all;

%% Parameters.
N = 300; % Input-Output characteristics / Binomial fitness distribution: number of trials.
gamma = 0.1937; % Rewiring probability.
q = 0.375; % Output characteristics / Binomial fitness distribution: success probability.
p = 0.075; % Input characteristics success probability.
load_eta = true;
phi = 0.6322; % Bonacich centrality parameter.
aggr_sic_level = 2;
drop_sic_missing = true;
dt = 0.0018; % Time increment.
T = 250; % Maximum number of iterations.
S = 100; % Number of realizations.
T_burn_in = 1; % Burn in time.
selected_year = 2015;

tic

rng(1,'twister')

%% Load firm id mappings data.
fid = fopen(['Data/firm_factset_ids_sic_nation.csv']);
dat = textscan(fid,repmat('%s ',1,4), 'delimiter', ';', 'headerlines', 1); % id; firm_factset_id; firm_sic; firm_nation
fclose(fid);
factset_firm_id = str2double(strtrim(dat{2}));
factset_firm_sic_str = strtrim(dat{3});
factset_firm_nation_str = strtrim(dat{4});
clear('dat')

%% Convert nation strings to integers.
factset_firm_nation = nan(length(factset_firm_nation_str),1);
factset_firm_nation_str_init_unique = unique(factset_firm_nation_str);
ind = find(strcmp(factset_firm_nation_str_init_unique,''));
factset_firm_nation_str_init_unique(ind) = [];
for i=1:length(factset_firm_nation_str_init_unique)
    if(~strcmp(factset_firm_nation_str_init_unique{i},''))
        ind = find(strcmp(factset_firm_nation_str_init_unique{i},factset_firm_nation_str));
        if(~isempty(ind))
            factset_firm_nation(ind) = i;
        end
    end
end
factset_firm_nation_init = factset_firm_nation;
factset_firm_nation_str_init = factset_firm_nation_str;
ind = find(strcmp(factset_firm_nation_str_init,''));
factset_firm_nation_str_init(ind) = [];

%% Aggregate SIC codes.
factset_firm_sic = nan(length(factset_firm_sic_str),1);
for i=1:length(factset_firm_sic_str)
    if(~strcmp(factset_firm_sic_str{i},'NaN'))
        dummy = factset_firm_sic_str{i};
        dummy = dummy(1:aggr_sic_level);
        if(~isempty(str2num(dummy)))
            factset_firm_sic(i) = str2num(dummy);
        end
    end
end

%% Load adjacency matrix.
disp(['Load adjacency matrix from file...'])
load (['./Data/' num2str(selected_year) '_adjacency_matrix.mat'])

%% Load firm fixed effects (fitness values).
if(load_eta)
    disp(['Load firm fixed effects (fitness values) from file...'])
    load (['./Data/mle_eta.mat']);
    eta = zeros(length(factset_firm_id),1);
    mle_eta_id = mle_eta(:,1);
    mle_eta_fit = mle_eta(:,2);
    for i=1:length(factset_firm_id)
        ind = find(factset_firm_id(i)==mle_eta_id);
        if(~isempty(ind))
            eta(i) = mle_eta_fit(ind(1));
        end
    end
    clear('mle_eta')
    mle_eta = [factset_firm_id eta];
else
    eta = zeros(length(factset_firm_id),1);
    for i=1:length(factset_firm_id)
        k = binornd(N,q);
        a = (1-p)^N;
        b = 1 + p/(q*(1-p));
        eta(i) = a*(b^k-1);
    end
    mle_eta = [factset_firm_id eta];
end
clear('eta')

%% Drop firms with missing SIC codes.
if(drop_sic_missing)
    disp(['Drop firms with missing SIC codes...'])
    ind = find(~isnan(factset_firm_sic));
    factset_firm_id = factset_firm_id(ind);
    factset_firm_sic = factset_firm_sic(ind);
    factset_firm_nation_str = factset_firm_nation_str(ind);
    factset_firm_nation = factset_firm_nation(ind);
    A = A(ind,ind);
    mle_eta = mle_eta(ind,:);
end

A_benchmark = A;

%% Estimated rewiring rates.
rewiring_gamma = 0.1937*ones(length(rewiring_sic),1);

%% Estimated exit rates.
lifetime_delta = 0.097479*ones(length(lifetime_sic),1); 
lifetime_zeta = 0.0060066*ones(length(lifetime_sic),1); 
lifetime_rho = 0.0082788*ones(length(lifetime_sic),1); 

%% Load empirical in-degree distribution.
disp(['Load empirical in-degree distribution from file...'])
load('./Data/in_deg_distribution.mat');
ind = find(in_deg_distribution(:,1)==selected_year);
emp_indeg_distribution = in_deg_distribution(end,2:end);

%% Load list of European countries.
disp('Load firms locations..')
fid = fopen(['./Data/European_countries.csv']); %
dat = textscan(fid,repmat('%q ',1,4), 'delimiter', ',', 'headerlines', 1); 
fclose(fid);
country_iso_2 = strtrim(dat{2});
country_eu_member = str2double(strtrim(dat{4}));
clear('dat')
ind = find(country_eu_member==0);
country_iso_2(ind) = []; % Keep only EU member states.
eu_iso_2 = country_iso_2;

%% Iterate over different scenarios.
disp(['Iterate over different scenarios.'])
scenario = {'benchmark'; 'autarky'};
n_scenario = [];
n_scenario_us = [];
n_scenario_asia = [];
n_scenario_eu = [];
bonacich_mean_scenario = [];
bonacich_mean_scenario_us = [];
bonacich_mean_scenario_asia = [];
bonacich_mean_scenario_eu = [];
clear('A')
n_time_save = [];
bonacich_time_save = [];
n_time_save_us = [];
bonacich_time_save_us = [];
n_time_save_asia = [];
bonacich_time_save_asia = [];
n_time_save_eu = [];
bonacich_time_save_eu = [];
for c=1:length(scenario)
    
    if(strcmp(scenario(c),'benchmark'))
        fprintf('Compute benchmark... \n');
    else
        fprintf('Compute counterfactual... \n');
    end
    A_0 = A_benchmark;
    n_0 = length(A_0);
    disp(['Number of firms ' num2str(n_0) '...']);
    
    %% Initialize vectors.
    n_time = zeros(T,S);
    n_time_us = zeros(T,S);
    n_time_asia = zeros(T,S);
    n_time_eu = zeros(T,S);
    bonacich_sum_time = zeros(T,S);
    bonacich_sum_time_us = zeros(T,S);
    bonacich_sum_time_asia = zeros(T,S);
    bonacich_sum_time_eu = zeros(T,S);
    
    %% Iterating over different realizations.
    for s = 1:S
        
        fprintf('Realization %3.0f.\n', s);
        
        %% Setting the seed for the random number generator.
        rand('seed',s);
        
        %% Generate initial network and assign initial fitness values.
        A = A_0;
        n = length(A);
        eta = mle_eta(:,2); % Assign initial fitness values.
        
        %% Assign delta, zeta and rho to incumbent firms.
        delta = 0.097479*ones(length(A_0),1); %0.0973*ones(length(A_0),1);
        zeta = 0.0060066*ones(length(A_0),1); %0.0037*ones(length(A_0),1);
        rho = 0.0082788*ones(length(A_0),1); %0.0065*ones(length(A_0),1);
        for i=1:length(A)
            ind = find(lifetime_sic==factset_firm_sic(i));
            if(~isempty(ind))
                delta(i) = lifetime_delta(ind);
                zeta(i) = lifetime_zeta(ind);
                rho(i) = lifetime_rho(ind);
            end
        end
        
        %% Assign the number of initial suppliers.
        num_init_suppl = full(sum(A,1)); % in-degree.
        
        %% Assign US indicator.
        is_us = zeros(length(A),1);
        ind_us = find(strcmp(factset_firm_nation_str,'US'));
        is_us(ind_us) = 1;
        
        %% Assign Asia indicator.
        is_china = zeros(length(A),1);
        ind_china = find(strcmp(factset_firm_nation_str,'CN'));
        is_china(ind_china) = 1;
        is_taiwan = zeros(length(A),1);
        ind_taiwan = find(strcmp(factset_firm_nation_str,'TW'));
        is_taiwan(ind_taiwan) = 1;
        is_japan = zeros(length(A),1);
        ind_japan = find(strcmp(factset_firm_nation_str,'JP'));
        is_japan(ind_japan) = 1;
        is_korea = zeros(length(A),1);
        ind_korea = find(strcmp(factset_firm_nation_str,'KR'));
        is_korea(ind_korea) = 1;
        is_singapore = zeros(length(A),1);
        ind_singapore = find(strcmp(factset_firm_nation_str,'SG'));
        is_singapore(ind_singapore) = 1;
        is_hong_kong = zeros(length(A),1);
        ind_hong_kong = find(strcmp(factset_firm_nation_str,'HK'));
        is_hong_kong(ind_hong_kong) = 1;
        is_asia = zeros(length(A),1);
        is_asia(ind_china) = 1;
        is_asia(ind_taiwan) = 1;
        is_asia(ind_japan) = 1;
        is_asia(ind_korea) = 1;
        is_asia(ind_singapore) = 1;
        is_asia(ind_hong_kong) = 1;
        ind_asia = [ind_china; ind_taiwan; ind_japan; ind_korea; ind_singapore; ind_hong_kong];
        firm_nation = factset_firm_nation;
        firm_nation_str = factset_firm_nation_str;
        
        %% Assign EU indicator.
        ind_eu = [];
        for j=1:length(eu_iso_2)
            ind_eu = [ind_eu; find(strcmp(factset_firm_nation_str,eu_iso_2{j}))];
        end
        is_eu = zeros(length(A),1);
        is_eu(ind_eu) = 1;
        
        %links_blocked = links_blocked_init;
        
        %% Network formation.
        for t=1:(T_burn_in+T)
            
            %if(mod(t,250)==0)
            %    fprintf('Network update %3.0f.\n', t);
            %end
            
            %% Removal all external links in the US, Europe and Eastern Asia.
            if(t==T_burn_in && strcmp(scenario(c),'autarky'))
                
                ind_us = find(is_us);
                ind_not_us = setdiff([1:1:length(A)],ind_us)'; % C = setdiff(A,B) returns the data in A that is not in B, with no repetitions. C is in sorted order.
                for k=1:length(ind_us)
                    A(ind_us(k),ind_not_us) = 0;
                    A(ind_not_us,ind_us(k)) = 0;
                end
                ind_eu = find(is_eu);
                ind_not_eu = setdiff([1:1:length(A)],ind_eu)';
                for k=1:length(ind_eu)
                    A(ind_eu(k),ind_not_eu) = 0;
                    A(ind_not_eu,ind_eu(k)) = 0;
                end
                ind_asia = find(is_asia);
                ind_not_asia = setdiff([1:1:length(A)],ind_asia)';
                for k=1:length(ind_asia)
                    A(ind_asia(k),ind_not_asia) = 0;
                    A(ind_not_asia,ind_asia(k)) = 0;
                end
            end
            
            %% Compute the time varying statistics.
            n = length(A);
            u = ones(n,1);
            out_d = full(sum(A,2)); % out-degree.
            in_d = full(sum(A,1))'; % in-degree.
            deg_av_time(t,s) = mean(in_d);
            in_d_std_time(t,s) = std(in_d);
            out_d_std_time(t,s) = std(out_d);
            n_time(t,s) = length(A);
            d_inv = 1./in_d;
            d_inv(isinf(d_inv)) = 1;
            W = A.*d_inv; % A.*d = A*diag(d); D_inv = diag(d); A./max(u,in_d); % A./d = A*inv(diag(d)).
            bonacich = u;
            %bonacich = inv(eye(n) - phi*W)*u;
            bonacich = (eye(n) - phi*W)\u; % Replace inv(A)*b with A\b.
            bonacich_sum_time(t,s) = nansum(bonacich);
            ind = find(is_us);
            n_time_us(t,s) = length(ind);
            bonacich_sum_time_us(t,s) = nansum(bonacich(ind));
            ind = find(is_asia);
            n_time_asia(t,s) = length(ind);
            bonacich_sum_time_asia(t,s) = nansum(bonacich(ind));
            ind = find(is_eu);
            n_time_eu(t,s) = length(ind);
            bonacich_sum_time_eu(t,s) = nansum(bonacich(ind));
            
            %% With probability delta*dt, remove each node at random.
            if(length(A)>1)
                ind_remove = [];
                for k=1:length(A)
                    if(rand<delta(k)*dt)
                        ind_remove = [ind_remove k];
                    end
                end
                while(~isempty(ind_remove))
                    k = ind_remove(1);
                    ind_remove(1) = [];
                    ind = find(ind_remove>k);
                    for j=1:length(ind)
                        ind_remove(ind(j)) = ind_remove(ind(j))-1;
                    end
                    A(k,:) = [];
                    A(:,k) = [];
                    eta(k) = [];
                    delta(k) = [];
                    zeta(k) = [];
                    rho(k) = [];
                    num_init_suppl(k) = [];
                    is_us(k) = [];
                    is_asia(k) = [];
                    is_eu(k) = [];
                    firm_nation(k) = [];
                    firm_nation_str(k) = [];
                end
            end
            
            %% With probability zeta*dt, remove each node at random if its after shock profits are negative (and profits are proportional to the out-degree).
            if(length(A)>1)
                ind_remove = [];
                for k=1:length(A)
                    if(rand<zeta(k)*dt/bonacich(k))
                        ind_remove = [ind_remove k];
                    end
                end
                while(~isempty(ind_remove))
                    k = ind_remove(1);
                    ind_remove(1) = [];
                    ind = find(ind_remove>k);
                    for j=1:length(ind)
                        ind_remove(ind(j)) = ind_remove(ind(j))-1;
                    end
                    A(k,:) = [];
                    A(:,k) = [];
                    eta(k) = [];
                    delta(k) = [];
                    zeta(k) = [];
                    rho(k) = [];
                    num_init_suppl(k) = [];
                    is_us(k) = [];
                    is_asia(k) = [];
                    is_eu(k) = [];
                    firm_nation(k) = [];
                    firm_nation_str(k) = [];
                end
            end
            
            %% With probability rho*dt select each firm which has an in-degree zero and remove it from the network.
            if(length(A)>1)
                ind_remove = [];
                for k=1:length(A)
                    if(in_d(k)==0 && rand<rho(k)*dt)
                        ind_remove = [ind_remove k];
                    end
                end
                while(~isempty(ind_remove))
                    k = ind_remove(1);
                    ind_remove(1) = [];
                    ind = find(ind_remove>k);
                    for j=1:length(ind)
                        ind_remove(ind(j)) = ind_remove(ind(j))-1;
                    end
                    A(k,:) = [];
                    A(:,k) = [];
                    eta(k) = [];
                    delta(k) = [];
                    zeta(k) = [];
                    rho(k) = [];
                    num_init_suppl(k) = [];
                    is_us(k) = [];
                    is_asia(k) = [];
                    is_eu(k) = [];
                    firm_nation(k) = [];
                    firm_nation_str(k) = [];
                    for l=1:length(ind_remove)
                        if(ind_remove(l)>k)
                            ind_remove(l) = ind_remove(l)-1;
                        end
                    end
                end
            end
            
            %% With probability gamma*dt select each firm which has in-degree less than the inital in-degree and find a new supplier from the incumbent firms (rewiring).
            if(length(A)>1)
                ind_rewire = [];
                for k=1:length(A)
                    if(in_d(k)< num_init_suppl(k) && rand<gamma*dt) %if(in_d(k)==0 && rand<gamma(k)*dt)
                        ind_rewire = [ind_rewire k];
                    end
                end
                n = length(A);
                while(~isempty(ind_rewire))
                    
                    ind = ind_rewire(1);
                    ind_rewire(1) = [];
                    
                    %% Fitness (eta) based attachment.
                    if(strcmp(scenario(c),'autarky') && t>T_burn_in)
                        if(is_us(ind))
                            dummy = find(is_us==1);
                        elseif(is_eu(ind))
                            dummy = find(is_eu==1);
                        elseif(is_asia(ind))
                            dummy = find(is_asia==1);
                        else
                            dummy = [find(is_us==1); find(is_eu==1); find(is_asia==1)];
                            dummy = setdiff([1:1:n],dummy)'; % C = setdiff(A,B) returns the data in A that is not in B, with no repetitions. C is in sorted order.
                        end
                        dummy(dummy==ind) = [];
                        j = 1;
                        f = eta(dummy)/length(dummy);
                        f = f/sum(f);
                        F = f(j);
                        u = rand;
                        while(F<u)
                            j = min(j+1,length(f));
                            F = F + f(j);
                        end
                        A(dummy(j),ind) = 1;
                    else
                        j = 1;
                        f = eta/n;
                        f = f/sum(f);
                        F = f(j);
                        u = rand;
                        while(F<u)
                            j = min(j+1,length(f));
                            F = F + f(j);
                        end
                        A(j,ind) = 1;
                    end
                    
                end
            end
            
            %% Add a new node.
            n = length(A);
            A = [A; sparse(1,n)]; % Expanding the adjacency matrix by one column and row of zeros.
            A = [A, sparse(n+1,1)];
            n = length(A);
            
            %% Assign new fitness value.
            k = binornd(N,q);
            a = (1-p)^N;
            b = 1 + p/(q*(1-p));
            eta(n) = a*(b^k-1);
            
            %% Assign rates.
            delta(n) = lifetime_delta(randi(length(lifetime_delta)));
            zeta(n) = lifetime_zeta(randi(length(lifetime_zeta)));
            rho(n) = lifetime_rho(randi(length(lifetime_rho)));
            
            %% Assign country.
            country = factset_firm_nation_str_init(randi(length(factset_firm_nation_str_init)));
            firm_nation_str(n) = country;
            if(strcmp(country,'US'))
                is_us = [is_us; 1];
                is_asia = [is_asia; 0];
                is_eu = [is_eu; 0];
            elseif(strcmp(country,'CN'))
                is_us = [is_us; 0];
                is_eu = [is_eu; 0];
                is_asia = [is_asia; 1];
            elseif(strcmp(country,'TW'))
                is_us = [is_us; 0];
                is_eu = [is_eu; 0];
                is_asia = [is_asia; 1];
            elseif(strcmp(country,'JP'))
                is_us = [is_us; 0];
                is_eu = [is_eu; 0];
                is_asia = [is_asia; 1];
            elseif(strcmp(country,'KR'))
                is_us = [is_us; 0];
                is_eu = [is_eu; 0];
                is_asia = [is_asia; 1];
            elseif(strcmp(country,'SG'))
                is_us = [is_us; 0];
                is_eu = [is_eu; 0];
                is_asia = [is_asia; 1];
            elseif(strcmp(country,'HK'))
                is_us = [is_us; 0];
                is_eu = [is_eu; 0];
                is_asia = [is_asia; 1];
            else
                found_eu = 0;
                for j=1:length(eu_iso_2)
                    if(strcmp(country,eu_iso_2{j}))
                        is_us = [is_us; 0];
                        is_asia = [is_asia; 0];
                        is_eu = [is_eu; 1];
                        found_eu = 1;
                    end
                end
                if(~found_eu)
                    is_us = [is_us; 0];
                    is_asia = [is_asia; 0];
                    is_eu = [is_eu; 0];
                end
            end
            ind = find(strcmp(country,factset_firm_nation_str_init_unique));
            if(~isempty(ind) && ~isempty(factset_firm_nation_str_init_unique{ind}))
                firm_nation(n) = ind;
            else
                firm_nation(n) = NaN;
            end
            
            %% Draw the number of suppliers of the firm.
            l = 1;
            F = emp_indeg_distribution(l);
            u = rand;
            while(F<u)
                l = l+1;
                F = F + emp_indeg_distribution(l);
            end
            num_init_suppl(n) = l; % Record the number of initial suppliers.
            
            %% Draw l suppliers from the incumbent firms.
            for k=1:num_init_suppl(n)
                
                %% Fitness (eta) based attachment.
                if(strcmp(scenario(c),'autarky') && t>T_burn_in)
                    if(is_us(n))
                        dummy = find(is_us==1);
                    elseif(is_eu(n))
                        dummy = find(is_eu==1);
                    elseif(is_asia(n))
                        dummy = find(is_asia==1);
                    else
                        dummy = [find(is_us==1); find(is_eu==1); find(is_asia==1)];
                        dummy = setdiff([1:1:n],dummy)'; % C = setdiff(A,B) returns the data in A that is not in B, with no repetitions. C is in sorted order.
                    end
                    dummy(dummy==n) = [];
                    j = 1;
                    f = eta(dummy)/length(dummy);
                    F = f(j);
                    u = rand;
                    while(F<u)
                        j = min(j+1,length(f));
                        F = F + f(j);
                    end
                    A(dummy(j),n) = 1;
                else
                    j = 1;
                    f = eta/n;
                    f = f/sum(f);
                    F = f(j);
                    u = rand;
                    while(F<u)
                        j = min(j+1,length(f));
                        F = F + f(j);
                    end
                    A(j,n) = 1;
                end
                
            end
            
        end
        
    end
    
    %% Save statistics for different scenarios.
    n_time_mean = nanmean(n_time(T_burn_in:end,:),2);
    n_scenario = [n_scenario, n_time_mean];
    n_time_mean_us = nanmean(n_time_us(T_burn_in:end,:),2);
    n_scenario_us = [n_scenario_us n_time_mean_us];
    n_time_mean_asia = nanmean(n_time_asia(T_burn_in:end,:),2);
    n_scenario_asia = [n_scenario_asia n_time_mean_asia];
    n_time_mean_eu = nanmean(n_time_eu(T_burn_in:end,:),2);
    n_scenario_eu = [n_scenario_eu n_time_mean_eu];
    
    bonacich_time_save = nanmean(bonacich_sum_time(T_burn_in:end,:),2);
    bonacich_mean_scenario = [bonacich_mean_scenario bonacich_time_save];
    bonacich_time_save_us = nanmean(bonacich_sum_time_us(T_burn_in:end,:),2);
    bonacich_mean_scenario_us = [bonacich_mean_scenario_us bonacich_time_save_us];
    bonacich_time_save_asia = nanmean(bonacich_sum_time_asia(T_burn_in:end,:),2);
    bonacich_mean_scenario_asia = [bonacich_mean_scenario_asia bonacich_time_save_asia];
    bonacich_time_save_eu = nanmean(bonacich_sum_time_eu(T_burn_in:end,:),2);
    bonacich_mean_scenario_eu = [bonacich_mean_scenario_eu bonacich_time_save_eu];
    
end

%% Plot the impulse response function for the US.
disp(['Plot the impulse response function for the US...']);
figure();
set(0, 'defaultTextInterpreter', 'latex');
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
benchmark = bonacich_mean_scenario_us(:,1);
counterfactual = bonacich_mean_scenario_us(:,2);
x = dt*[1:1:length(benchmark)]';
y1 = benchmark;
y2 = counterfactual;
ind = find(y1>0);
x = x(ind);
y1 = y1(ind);
y2 = y2(ind);
y2 = 100*(y2-benchmark(1))./benchmark(1);
%y2 = 100*(y2-counterfactual(1))./counterfactual(1);
%plot(x,y2,'-r','LineWidth', 2);
hold on
y1 = 100*(y1-benchmark(1))./benchmark(1);
%plot(x,y1,'-k','LineWidth', 2);
hold on
h = fill([x;flipud(x)],[y1;flipud(y2)],'r','LineWidth',0.1);
set(h,'facealpha',0.25)
hold on
plot(x,y1,'-k','LineWidth', 2); % benchmark
hold on
plot(x,y2,'-r','LineWidth', 2); % counterfactual with autarky.
hold on
axis tight
hold on
%set(gca, 'yscale', 'log')
%set(gca, 'xscale', 'log')
xlabel('$t$');
ylabel('$(Y_t-Y_0) / Y_0$ \%');
title('US Firms')
set(gca,'FontSize',20);
filename = strcat('Figures/','Prodnet_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(nanmean(delta)),'_zeta_',num2str(nanmean(zeta)),'_gamma_',num2str(nanmean(gamma)),'_rho_',num2str(nanmean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_impulse_response_autarky_us_eu_east_asia_us.eps');
print('-depsc',filename)

%% Plot the impulse response function for the Asia.
disp(['Plot the impulse response function for the Asia...']);
figure();
set(0, 'defaultTextInterpreter', 'latex');
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
benchmark = bonacich_mean_scenario_asia(:,1);
counterfactual = bonacich_mean_scenario_asia(:,2);
x = dt*[1:1:length(benchmark)]';
y1 = benchmark;
y2 = counterfactual;
ind = find(y1>0);
x = x(ind);
y1 = y1(ind);
y2 = y2(ind);
y2 = 100*(y2-benchmark(1))./benchmark(1);
%y2 = 100*(y2-counterfactual(1))./counterfactual(1);
%plot(x,y2,'-r','LineWidth', 2);
hold on
y1 = 100*(y1-benchmark(1))./benchmark(1);
%plot(x,y1,'-k','LineWidth', 2);
hold on
h = fill([x;flipud(x)],[y1;flipud(y2)],'r','LineWidth',0.1);
set(h,'facealpha',0.25)
hold on
plot(x,y1,'-k','LineWidth', 2); % benchmark
hold on
plot(x,y2,'-r','LineWidth', 2); % counterfactual with autarky.
hold on
axis tight
%set(gca, 'yscale', 'log')
%set(gca, 'xscale', 'log')
xlabel('$t$');
ylabel('$(Y_t-Y_0) / Y_0$ \%');
title('Asian Firms')
set(gca,'FontSize',20);
filename = strcat('Figures/','Prodnet_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(nanmean(delta)),'_zeta_',num2str(nanmean(zeta)),'_gamma_',num2str(nanmean(gamma)),'_rho_',num2str(nanmean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_impulse_response_autarky_us_eu_east_asia_asia.eps');
print('-depsc',filename)

%% Plot the impulse response function for the EU.
disp(['Plot the impulse response function for the EU...']);
figure();
set(0, 'defaultTextInterpreter', 'latex');
set(gca,'Layer','top')
set(gca,'FontSize',20);
box on
benchmark = bonacich_mean_scenario_eu(:,1);
counterfactual = bonacich_mean_scenario_eu(:,2);
x = dt*[1:1:length(benchmark)]';
y1 = benchmark;
y2 = counterfactual;
ind = find(y1>0);
x = x(ind);
y1 = y1(ind);
y2 = y2(ind);
y2 = 100*(y2-benchmark(1))./benchmark(1);
%y2 = 100*(y2-counterfactual(1))./counterfactual(1);
%plot(x,y2,'-r','LineWidth', 2);
hold on
y1 = 100*(y1-benchmark(1))./benchmark(1);
%plot(x,y1,'-k','LineWidth', 2);
hold on
h = fill([x;flipud(x)],[y1;flipud(y2)],'r','LineWidth',0.1);
set(h,'facealpha',0.25)
hold on
plot(x,y1,'-k','LineWidth', 2); % benchmark
hold on
plot(x,y2,'-r','LineWidth', 2); % counterfactual with autarky.
hold on
axis tight
%set(gca, 'yscale', 'log')
%set(gca, 'xscale', 'log')
xlabel('$t$');
ylabel('$(Y_t-Y_0) / Y_0$ \%');
title('EU Firms')
set(gca,'FontSize',20);
filename = strcat('Figures/','Prodnet_N_',num2str(N),'_p_',num2str(p),'_q_',num2str(q),'_delta_',num2str(nanmean(delta)),'_zeta_',num2str(nanmean(zeta)),'_gamma_',num2str(nanmean(gamma)),'_rho_',num2str(nanmean(rho)),'_T_',num2str(T),'_S_',num2str(S),'_impulse_response_autarky_us_eu_east_asia_eu.eps');
print('-depsc',filename)

t = toc;
disp(['Elapsed time is ' datestr(datenum(0,0,0,0,0,t),'HH:MM:SS')])